home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 60 / IOPROG_60.ISO / soft / c++ / gsl-1.1.1-setup.exe / {app} / src / eigen / test.c < prev   
Encoding:
C/C++ Source or Header  |  2002-04-18  |  11.4 KB  |  390 lines
  1. /* eigen/test.c
  2.  * 
  3.  * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
  4.  * 
  5.  * This program is free software; you can redistribute it and/or modify
  6.  * it under the terms of the GNU General Public License as published by
  7.  * the Free Software Foundation; either version 2 of the License, or (at
  8.  * your option) any later version.
  9.  * 
  10.  * This program is distributed in the hope that it will be useful, but
  11.  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13.  * General Public License for more details.
  14.  * 
  15.  * You should have received a copy of the GNU General Public License
  16.  * along with this program; if not, write to the Free Software
  17.  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  */
  19.  
  20. /* Author:  G. Jungman
  21.  */
  22.  
  23. #include <config.h>
  24. #include <stdlib.h>
  25. #include <gsl/gsl_test.h>
  26. #include <gsl/gsl_math.h>
  27. #include <gsl/gsl_blas.h>
  28. #include <gsl/gsl_ieee_utils.h>
  29. #include <gsl/gsl_complex_math.h>
  30. #include <gsl/gsl_eigen.h>
  31.  
  32. gsl_matrix *
  33. create_hilbert_matrix(int size)
  34. {
  35.   int i, j;
  36.   gsl_matrix * m = gsl_matrix_alloc(size, size);
  37.   for(i=0; i<size; i++) {
  38.     for(j=0; j<size; j++) {
  39.       gsl_matrix_set(m, i, j, 1.0/(i+j+1.0));
  40.     }
  41.   }
  42.   return m;
  43. }
  44.  
  45. gsl_matrix *
  46. create_random_symm_matrix(int size)
  47. {
  48.   int i, j;
  49.   unsigned long k = 1;
  50.   gsl_matrix * m = gsl_matrix_alloc(size, size);
  51.   for(i=0; i<size; i++) {
  52.     for(j=i; j<size; j++) {
  53.       double x;
  54.       k = (69069 * k + 1) & 0xffffffffUL;
  55.       x = k / 4294967296.0;
  56.       gsl_matrix_set(m, i, j, x);
  57.       gsl_matrix_set(m, j, i, x);
  58.     }
  59.   }
  60.   return m;
  61. }
  62.  
  63. gsl_matrix_complex *
  64. create_random_herm_matrix(int size)
  65. {
  66.   int i, j;
  67.   unsigned long k = 1;
  68.   gsl_matrix_complex * m = gsl_matrix_complex_alloc(size, size);
  69.   for(i=0; i<size; i++) {
  70.     for(j=i; j<size; j++) {
  71.       gsl_complex z;
  72.       k = (69069 * k + 1) & 0xffffffffUL;
  73.       GSL_REAL(z) = k / 4294967296.0;
  74.       k = (69069 * k + 1) & 0xffffffffUL;
  75.       GSL_IMAG(z) = (i == j) ? 0 : k / 4294967296.0;
  76.       gsl_matrix_complex_set(m, i, j, z);
  77.       gsl_matrix_complex_set(m, j, i, gsl_complex_conjugate(z));
  78.     }
  79.   }
  80.   return m;
  81. }
  82.  
  83. void
  84. test_eigen_results (size_t N, const gsl_matrix * m, const gsl_vector * eval, 
  85.                     const gsl_matrix * evec, const char * desc,
  86.                     const char * desc2)
  87. {
  88.   size_t i,j;
  89.  
  90.   gsl_vector * x = gsl_vector_alloc(N);
  91.   gsl_vector * y = gsl_vector_alloc(N);
  92.  
  93.   /* check eigenvalues */
  94.  
  95.   for (i = 0; i < N; i++)
  96.     {
  97.       double ei = gsl_vector_get (eval, i);
  98.       gsl_vector_const_view vi = gsl_matrix_const_column(evec, i);
  99.       gsl_vector_memcpy(x, &vi.vector);
  100.       /* compute y = m x (should = lambda v) */
  101.       gsl_blas_dgemv (CblasNoTrans, 1.0, m, x, 0.0, y);
  102.       for (j = 0; j < N; j++)
  103.         {
  104.           double xj = gsl_vector_get (x, j);
  105.           double yj = gsl_vector_get (y, j);
  106.           gsl_test_rel(yj, ei * xj,  1e8 * GSL_DBL_EPSILON, 
  107.                        "%s, eigenvalue(%d,%d), %s", desc, i, j, desc2);
  108.         }
  109.     }
  110.  
  111.   /* check eigenvectors are orthonormal */
  112.  
  113.   for (i = 0; i < N; i++)
  114.     {
  115.       gsl_vector_const_view vi = gsl_matrix_const_column(evec, i);
  116.       double nrm_v = gsl_blas_dnrm2(&vi.vector);
  117.       gsl_test_rel (nrm_v, 1.0, N * GSL_DBL_EPSILON, "%s, normalized(%d), %s", 
  118.                     desc, i, desc2);
  119.     }
  120.  
  121.   for (i = 0; i < N; i++)
  122.     {
  123.       gsl_vector_const_view vi = gsl_matrix_const_column(evec, i);
  124.       for (j = i + 1; j < N; j++)
  125.         {
  126.           gsl_vector_const_view vj = gsl_matrix_const_column(evec, j);
  127.           double vivj;
  128.           gsl_blas_ddot (&vi.vector, &vj.vector, &vivj);
  129.           gsl_test_abs (vivj, 0.0, N * GSL_DBL_EPSILON, 
  130.                         "%s, orthogonal(%d,%d), %s", desc, i, j, desc2);
  131.         }
  132.     }
  133.  
  134.   gsl_vector_free(x);
  135.   gsl_vector_free(y);
  136. }
  137.  
  138.  
  139. void
  140. test_eigenvalues (size_t N, const gsl_vector *eval, const gsl_vector * eval2, 
  141.                   const char * desc, const char * desc2)
  142. {
  143.   size_t i;
  144.   for (i = 0; i < N; i++)
  145.     {
  146.       double ei = gsl_vector_get (eval, i);
  147.       double e2i = gsl_vector_get (eval2, i);
  148.       gsl_test_rel(ei, e2i, GSL_DBL_EPSILON, "%s, direct eigenvalue(%d), %s",
  149.                    desc, i, desc2);
  150.     }
  151. }
  152.  
  153.  
  154. void
  155. test_eigen_complex_results (size_t N, const gsl_matrix_complex * m, 
  156.                             const gsl_vector * eval, 
  157.                             const gsl_matrix_complex * evec, 
  158.                             const char * desc,
  159.                             const char * desc2)
  160. {
  161.   size_t i,j;
  162.  
  163.   gsl_vector_complex * x = gsl_vector_complex_alloc(N);
  164.   gsl_vector_complex * y = gsl_vector_complex_alloc(N);
  165.  
  166.   /* check eigenvalues */
  167.  
  168.   for (i = 0; i < N; i++)
  169.     {
  170.       double ei = gsl_vector_get (eval, i);
  171.       gsl_vector_complex_const_view vi = gsl_matrix_complex_const_column(evec, i);
  172.       gsl_vector_complex_memcpy(x, &vi.vector);
  173.       /* compute y = m x (should = lambda v) */
  174.       gsl_blas_zgemv (CblasNoTrans, GSL_COMPLEX_ONE, m, x, 
  175.                       GSL_COMPLEX_ZERO, y);
  176.       for (j = 0; j < N; j++)
  177.         {
  178.           gsl_complex xj = gsl_vector_complex_get (x, j);
  179.           gsl_complex yj = gsl_vector_complex_get (y, j);
  180.           gsl_test_rel(GSL_REAL(yj), ei * GSL_REAL(xj), 1e8*GSL_DBL_EPSILON, 
  181.                        "%s, eigenvalue(%d,%d), real, %s", desc, i, j, desc2);
  182.           gsl_test_rel(GSL_IMAG(yj), ei * GSL_IMAG(xj), 1e8*GSL_DBL_EPSILON, 
  183.                        "%s, eigenvalue(%d,%d), imag, %s", desc, i, j, desc2);
  184.         }
  185.     }
  186.  
  187.   /* check eigenvectors are orthonormal */
  188.  
  189.   for (i = 0; i < N; i++)
  190.     {
  191.       gsl_vector_complex_const_view vi = gsl_matrix_complex_const_column(evec, i);
  192.       double nrm_v = gsl_blas_dznrm2(&vi.vector);
  193.       gsl_test_rel (nrm_v, 1.0, N * GSL_DBL_EPSILON, "%s, normalized(%d), %s", 
  194.                     desc, i, desc2);
  195.     }
  196.  
  197.   for (i = 0; i < N; i++)
  198.     {
  199.       gsl_vector_complex_const_view vi = gsl_matrix_complex_const_column(evec, i);
  200.       for (j = i + 1; j < N; j++)
  201.         {
  202.           gsl_vector_complex_const_view vj 
  203.             = gsl_matrix_complex_const_column(evec, j);
  204.           gsl_complex vivj;
  205.           gsl_blas_zdotc (&vi.vector, &vj.vector, &vivj);
  206.           gsl_test_abs (gsl_complex_abs(vivj), 0.0, N * GSL_DBL_EPSILON, 
  207.                         "%s, orthogonal(%d,%d), %s", desc, i, j, desc2);
  208.         }
  209.     }
  210.  
  211.   gsl_vector_complex_free(x);
  212.   gsl_vector_complex_free(y);
  213. }
  214.  
  215.  
  216. void
  217. test_eigen_symm(const char * desc, const gsl_matrix * m)
  218. {
  219.   size_t N = m->size1;
  220.  
  221.   gsl_matrix * A = gsl_matrix_alloc(N, N);
  222.   gsl_matrix * evec = gsl_matrix_alloc(N, N);
  223.   gsl_vector * eval = gsl_vector_alloc(N);
  224.   gsl_vector * eval2 = gsl_vector_alloc(N);
  225.  
  226.   gsl_eigen_symm_workspace * w1 = gsl_eigen_symm_alloc (N);
  227.   gsl_eigen_symmv_workspace * w2 = gsl_eigen_symmv_alloc (N);
  228.  
  229.   gsl_matrix_memcpy(A, m);
  230.   gsl_eigen_symmv(A, eval, evec, w2);
  231.   test_eigen_results (N, m, eval, evec, desc, "unsorted");
  232.  
  233.   gsl_matrix_memcpy(A, m);
  234.   gsl_eigen_symm(A, eval2, w1);
  235.   test_eigenvalues (N, eval, eval2, desc, "unsorted");
  236.  
  237.   gsl_eigen_symmv_sort (eval, evec, GSL_EIGEN_SORT_VAL_ASC);
  238.   test_eigen_results (N, m, eval, evec, desc, "val/asc");
  239.  
  240.   gsl_eigen_symmv_sort (eval, evec, GSL_EIGEN_SORT_VAL_DESC);
  241.   test_eigen_results (N, m, eval, evec, desc, "val/desc");
  242.  
  243.   gsl_eigen_symmv_sort (eval, evec, GSL_EIGEN_SORT_ABS_ASC);
  244.   test_eigen_results (N, m, eval, evec, desc, "abs/asc");
  245.  
  246.   gsl_eigen_symmv_sort (eval, evec, GSL_EIGEN_SORT_VAL_DESC);
  247.   test_eigen_results (N, m, eval, evec, desc, "abs/desc");
  248.  
  249.   gsl_eigen_symm_free (w1);
  250.   gsl_eigen_symmv_free (w2);
  251.  
  252.   gsl_matrix_free(A);
  253.   gsl_matrix_free(evec);
  254.   gsl_vector_free(eval);
  255.   gsl_vector_free(eval2);
  256. }
  257.  
  258.  
  259. void
  260. test_eigen_herm(const char * desc, const gsl_matrix_complex * m)
  261. {
  262.   size_t N = m->size1;
  263.  
  264.   gsl_matrix_complex * A = gsl_matrix_complex_alloc(N, N);
  265.   gsl_matrix_complex * evec = gsl_matrix_complex_alloc(N, N);
  266.   gsl_vector * eval = gsl_vector_alloc(N);
  267.   gsl_vector * eval2 = gsl_vector_alloc(N);
  268.  
  269.   gsl_eigen_herm_workspace * w1 = gsl_eigen_herm_alloc (N);
  270.   gsl_eigen_hermv_workspace * w2 = gsl_eigen_hermv_alloc (N);
  271.  
  272.   gsl_matrix_complex_memcpy(A, m);
  273.   gsl_eigen_hermv(A, eval, evec, w2);
  274.   test_eigen_complex_results (N, m, eval, evec, desc, "unsorted");
  275.  
  276.   gsl_matrix_complex_memcpy(A, m);
  277.   gsl_eigen_herm(A, eval2, w1);
  278.   test_eigenvalues (N, eval, eval2, desc, "unsorted");
  279.  
  280.   gsl_eigen_hermv_sort (eval, evec, GSL_EIGEN_SORT_VAL_ASC);
  281.   test_eigen_complex_results (N, m, eval, evec, desc, "val/asc");
  282.  
  283.   gsl_eigen_hermv_sort (eval, evec, GSL_EIGEN_SORT_VAL_DESC);
  284.   test_eigen_complex_results (N, m, eval, evec, desc, "val/desc");
  285.  
  286.   gsl_eigen_hermv_sort (eval, evec, GSL_EIGEN_SORT_ABS_ASC);
  287.   test_eigen_complex_results (N, m, eval, evec, desc, "abs/asc");
  288.  
  289.   gsl_eigen_hermv_sort (eval, evec, GSL_EIGEN_SORT_VAL_DESC);
  290.   test_eigen_complex_results (N, m, eval, evec, desc, "abs/desc");
  291.  
  292.   gsl_eigen_herm_free (w1);
  293.   gsl_eigen_hermv_free (w2);
  294.  
  295.   gsl_matrix_complex_free(A);
  296.   gsl_matrix_complex_free(evec);
  297.   gsl_vector_free(eval);
  298.   gsl_vector_free(eval2);
  299. }
  300.  
  301.  
  302. void
  303. test_eigen_jacobi(const char * desc, const gsl_matrix * m)
  304. {
  305.   size_t N = m->size1;
  306.   unsigned int nrot;
  307.  
  308.   gsl_matrix * A = gsl_matrix_alloc(N, N);
  309.   gsl_matrix * evec = gsl_matrix_alloc(N, N);
  310.   gsl_vector * eval = gsl_vector_alloc(N);
  311.  
  312.   gsl_matrix_memcpy(A, m);
  313.   gsl_eigen_jacobi(A, eval, evec, 1000, &nrot);
  314.   gsl_eigen_symmv_sort(eval, evec, GSL_EIGEN_SORT_VAL_ASC);
  315.  
  316.   test_eigen_results (N, m, eval, evec, desc, "");
  317.  
  318.   gsl_matrix_free(A);
  319.   gsl_matrix_free(evec);
  320.   gsl_vector_free(eval);
  321. }
  322.  
  323.  
  324. int test_invert_jacobi(void)
  325. {
  326.   int s = 0;
  327.   int i, j;
  328.   gsl_matrix * hminv = gsl_matrix_alloc(10, 10);
  329.   gsl_matrix * id    = gsl_matrix_alloc(10, 10);
  330.  
  331.   /* 10x10 Hilbert matrix */
  332.   gsl_matrix * hm = create_hilbert_matrix(10);
  333.   gsl_eigen_invert_jacobi(hm, hminv, 1000);
  334.  
  335.   /* gsl_linalg_matmult(hm, hminv, id); */
  336.   gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, 1.0, hm, hminv, 0.0, id);
  337.  
  338.   for(i=0; i<10; i++) {
  339.     for(j=0; j<10; j++) {
  340.       double delta_ij = ( i == j ? 1.0 : 0.0 );
  341.       double id_ij    = gsl_matrix_get(id, i, j);
  342.       int rs = ( fabs(id_ij - delta_ij) > 5.0e-3 );
  343.       s += rs;
  344.     }
  345.   }
  346.  
  347.   gsl_matrix_free(hm);
  348.   gsl_matrix_free(hminv);
  349.   gsl_matrix_free(id);
  350.  
  351.   return s;
  352. }
  353.  
  354. int main()
  355. {
  356.   gsl_matrix *rs10 = create_random_symm_matrix (10);
  357.   gsl_matrix_complex *rh10 = create_random_herm_matrix (10);
  358.  
  359.   double r[] =  { 0,  0, -1,  0, 
  360.                   0,  1,  0,  1,
  361.                  -1,  0,  0,  0,
  362.                   0,  1,  0,  0 };
  363.  
  364.   double c[] =  { 0,0,  0,0, -1,0,  0,0, 
  365.                   0,0,  1,0,  0,0,  1,0,
  366.                  -1,0,  0,0,  0,0,  0,0,
  367.                   0,0,  1,0,  0,0,  0,0 };
  368.  
  369.   gsl_matrix_view s4 = gsl_matrix_view_array (r, 4, 4);
  370.   gsl_matrix_complex_view h4 = gsl_matrix_complex_view_array (c, 4, 4);
  371.  
  372.   gsl_ieee_env_setup ();
  373.  
  374.   test_eigen_symm("symm(4)", &s4.matrix);
  375.   test_eigen_herm("herm(4)", &h4.matrix);
  376.  
  377.   test_eigen_symm("symm(10)", rs10);
  378.   test_eigen_herm("herm(10)", rh10);
  379.  
  380.   /* gsl_matrix *h5 = create_hilbert_matrix (5); */
  381.   /* test_eigen_jacobi("hilbert(5)", h5); */
  382.   /* test_invert_jacobi(); */
  383.   /* gsl_matrix_free (h5); */
  384.  
  385.   gsl_matrix_complex_free (rh10);
  386.   gsl_matrix_free (rs10);
  387.  
  388.   exit (gsl_test_summary());
  389. }
  390.